clear
set matsize 11000
cd <directory>
**********************************************************************************************************************
** 													Figures
**********************************************************************************************************************


use "netatmo_u10and12and23_until10aug2015.dta", clear
g year = year(date)
g month = month(date)
g dow = dow(date)
g outdoorsq = OutdoorTemperature^2

set more off
foreach var of varlist Temperature OutdoorTemperature {
reg `var'  i.unit#i.year i.unit#i.month i.dow
predict `var'res, res
}

* FIGURE: INDOOR OUTDOOR GRADIENT
preserve
keep if OutdoorTemperature >=21.2 & OutdoorTemperature<=32.6
collapse (mean) Temperature (sem) sem=Temperature, by(OutdoorTemperature)
g cihigh = Temperature + 1.96*sem
g cilow = Temperature - 1.96*sem
twoway (scatter Temperature OutdoorTemperature, mcolor(blue) msize(tiny)) ///
(lpoly Temperature OutdoorTemperature, clcolor(black)) ///
(lpoly cihigh OutdoorTemperature, clpattern(dot) clcolor(black)) ///
(lpoly cilow OutdoorTemperature, clpattern(dot) clcolor(black)), ///
xline(24.9) xline(28.5)  ///
legend(order(2 "Locally-weighted polynomial smoothed gradient" 3 " 95% CIs")) ///
 ytitle("Indoor Temperature") xtitle("Outdoor Temperature") ///
note("Outdoor temperature trimmed at 1st and 99th percentiles. Locally-weighted polynomial smoothing used in regression fit.", size(vsmall)) ///
caption("Scatter depicts mean indoor temperature by outdoor temperature bins of .1 width. Vertical lines depict 25th and 75th percentiles.", size(vsmall))
	graph export "IndoorOutdoor.pdf", replace
	
	twoway (scatter Temperature OutdoorTemperature, mcolor(black) msize(tiny)) ///
(lpoly Temperature OutdoorTemperature, clcolor(black)) ///
(lpoly cihigh OutdoorTemperature, clpattern(dot) clcolor(black)) ///
(lpoly cilow OutdoorTemperature, clpattern(dot) clcolor(black)), ///
xline(24.9) xline(28.5)  ///
legend(order(2 "Locally-weighted polynomial smoothed gradient" 3 " 95% CIs")) ///
 ytitle("Indoor Temperature") xtitle("Outdoor Temperature") ///
note("Outdoor temperature trimmed at 1st and 99th percentiles. Locally-weighted polynomial smoothing used in regression fit.", size(vsmall)) ///
caption("Scatter depicts mean indoor temperature by outdoor temperature bins of .1 width. Vertical lines depict 25th and 75th percentiles.", size(vsmall))  scheme(s2mono)

		graph export "IndoorOutdoor_bw.pdf", replace
restore

* FIGURE 2A: Productivity temp gradient (pre-led)- raw figures

use "LEDprojecwithtamudatanew.dta", clear
preserve
keep if led==0
sum wbgtmeant, d
keep if wbgtmeant >=r(p1) & wbgtmeant<=r(p99)
replace wbgtmeant = round(wbgtmeant,0.1)
bys wbgtmeant led: egen N=count(mactualeff)
bys led: sum N
egen sem=sd(mactualeff)
replace sem=sem/sqrt(N)
g cihigh = mactualeff + 1.96*sem
g cilow = mactualeff - 1.96*sem
bys wbgtmeant led: egen scattermean=mean(mactualeff)
twoway (lpoly mactualeff wbgtmeant if led==0, clcolor(black) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(scatter scattermean wbgtmeant if led==0, mcolor(blue) msize(tiny) yscale(range(40(5)60))) ///
, scheme(s2color) legend(order(1 "Locally-weighted polynomial smoothed gradient" 2 " 95% CIs")) ///
xline(19)  ///
xtitle("Wet Bulb Globe Temperature") ytitle("Actual Efficiency", height(7)) ///
note("Scatter depicts mean efficiency residual by temperature residual bins of .1 width. Temperature residuals are trimmed at the ", size(vsmall)) ///
caption("1st and 99th. Vertical line depicts spline node.", size(vsmall))
	graph export "efftempgradient_noled_wbgt_RAW.pdf", replace
	
twoway (lpoly mactualeff wbgtmeant if led==0, clcolor(black) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(scatter scattermean wbgtmeant if led==0, mcolor(black) msize(tiny) yscale(range(40(5)60))) ///
,  scheme(s2mono) legend(order(1 "Locally-weighted polynomial smoothed gradient" 2 " 95% CIs")) ///
xline(19)  ///
xtitle("Wet Bulb Globe Temperature") ytitle("Actual Efficiency", height(7)) ///
note("Scatter depicts mean efficiency residual by temperature residual bins of .1 width. Temperature residuals are trimmed at the ", size(vsmall)) ///
caption("1st and 99th. Vertical line depicts spline node.", size(vsmall))

	graph export "efftempgradient_noled_wbgt_RAW_bw.pdf", replace
restore


* FIGURE: efficiency temp gradient by led - raw data
preserve
sum wbgtmeant, d
keep if wbgtmeant >=r(p1) & wbgtmeant<=r(p99)
replace wbgtmeant = round(wbgtmeant,0.1)
bys wbgtmeant led: egen N=count(mactualeff)
bys led: sum N
egen sem=sd(mactualeff)
replace sem=sem/sqrt(N)
g cihigh = mactualeff + 1.3723*sem
g cilow = mactualeff - 1.3723*sem
twoway (lpoly mactualeff wbgtmeant if led==0, clpattern(dash) clcolor(red) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==0, clpattern(dot) clcolor(red) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==0, clpattern(dot) clcolor(red) degree(1) bw(1)) ///
(lpoly mactualeff wbgtmeant if led==1, clcolor(blue) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==1, clpattern(dot) clcolor(blue) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==1, clpattern(dot) clcolor(blue) degree(1) bw(1)) ///
, scheme(s2color) legend(order(1 "Before LED" 4 "After LED")) ///
xline(19) xtitle("Wet Bulb Globe Temperature") ytitle("Actual Efficiency", height(7)) ///
note("Temperature trimmed at the 1st and 99th. Vertical line depicts spline node.", size(vsmall))
	graph export "efftempgradient_byled_wbgt_RAW.pdf", replace	
	
	
twoway (lpoly mactualeff wbgtmeant if led==0, clpattern(dash) clcolor(black) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==0, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(lpoly mactualeff wbgtmeant if led==1, clcolor(black) degree(1) bw(1)) ///
(lpoly cihigh wbgtmeant if led==1, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
(lpoly cilow wbgtmeant if led==1, clpattern(dot) clcolor(black) degree(1) bw(1)) ///
, scheme(s2mono) legend(order(1 "Before LED" 4 "After LED")) ///
xline(19) xtitle("Wet Bulb Globe Temperature") ytitle("Actual Efficiency", height(7)) ///
note("Temperature trimmed at the 1st and 99th. Vertical line depicts spline node.", size(vsmall))	
graph export "efftempgradient_byled_wbgt_RAW_bw.pdf", replace
restore



******************************************Appendix Figures***********************************************
use "weather_2009to2014",clear
collapse (mean) maxtemperature mintemperature wind relativehumidity precipitation solar wbgtmeant himeant, by(latitude longitude year month)
tab year
tempfile temp
save `temp', replace
* FIGURE:Cumulative Proportion of Factories Adopting LED 

use "LEDprojecwithtamudatanew.dta" , clear
keep unit yearled monthled latitude longitude
duplicates drop unit longitude latitude, force
ren monthled month
ren yearled year
drop if year==.
merge m:1 longitude latitude year month using `temp'
drop if _merge==2
cou
g n=1
collapse (sum) n (mean) wbgtmeant, by(year month)
g prop=n/30
gen timem = ym(year, month)

sort year month
g cumulativeprop=sum(prop)
format timem $dofm
label var cumulativeprop "Cumulative Share Adopting"
label var wbgtmeant "Mean WBGT"
label var timem "Month of LED Adoption"

twoway line cumulativeprop timem, yaxis(1) ytitle("Cumulative Proportion of Factories Adopting LED")|| line wbgtmeant  timem, lpattern(dash) yaxis(2) ytitle("Mean Temperature (WBGT)", axis(2))  xlabel(597(6)636)
graph export "Pictures\cumulativeprop_wbgtmeant.pdf", as(pdf) replace


twoway line cumulativeprop timem,   yaxis(1) ytitle("Cumulative Proportion of Factories Adopting LED")|| line wbgtmeant  timem, yaxis(2) ytitle("Mean Temperature (WBGT)", axis(2))  xlabel(597(6)636) scheme(s2mono)
graph export "Pictures\cumulativeprop_wbgtmeant_bw.pdf", as(pdf) replace

* FIGURE: Semiparametric differences and t-stats

use "LEDprojecwithtamudatanew.dta", clear

preserve 
sum wbgtmeantres, d
keep if wbgtmeantres >=r(p1) & wbgtmeantres<=r(p99)
replace wbgtmeantres = round(wbgtmeantres,0.1)
g ledbin = ledres>=0
g obs = _N
collapse mactualeffres probwbgtmeantres obs, by(wbgtmeantres led)
lpoly mactualeffres wbgtmeantres if led==0, generate(smooth0) at(wbgtmeantres) se(se0)
lpoly mactualeffres wbgtmeantres if led==1, generate(smooth1)  at(wbgtmeantres) se(se1)
g diff = smooth1-smooth0
g sem = sqrt(se1^2+se0^2)
g t=diff/sem
g cihigh = diff+1.96*sem
g cilow = diff-1.96*sem
twoway (connected diff wbgtmeantres, yaxis(1)) ///
(line cihigh wbgtmeantres, clcolor(black) clpattern(dash) yaxis(1)) ///
(line cilow wbgtmeantres, clcolor(black) clpattern(dash) yaxis(1)) ///
(lowess probwbgtmeantres wbgtmeantres, clpattern(dot) yaxis(2)), yline(0) ///
legend(order(1 "Effect Estimate" 2 "95% CIs" 4 "Temperature Density")) ///
ytitle("After LED - Before LED", axis(1)) ytitle("Temperature Density", axis(2)) xtitle("Outdoor Wet Bulb Globe Residuals") ///
caption("Horizontal line depicts difference of 0. Temperature is trimmed at the 1st and 99th percentiles.", size(vsmall))
	graph export "eff_semipardiff_wbgt.pdf", replace	
	
	twoway (connected diff wbgtmeantres, yaxis(1)) ///
(line cihigh wbgtmeantres, clcolor(black) clpattern(dash) yaxis(1)) ///
(line cilow wbgtmeantres, clcolor(black) clpattern(dash) yaxis(1)) ///
(lowess probwbgtmeantres wbgtmeantres, clpattern(dot) yaxis(2)), yline(0) ///
legend(order(1 "Effect Estimate" 2 "95% CIs" 4 "Temperature Density")) ///
ytitle("After LED - Before LED", axis(1)) ytitle("Temperature Density", axis(2)) xtitle("Outdoor Wet Bulb Globe Residuals") ///
caption("Horizontal line depicts difference of 0. Temperature is trimmed at the 1st and 99th percentiles.", size(vsmall)) scheme(s2mono)	
	graph export "eff_semipardiff_wbgt_bw.pdf", replace 



* FIGURE: Temperature distribution by led
preserve
sum wbgtmeantres, d
keep if wbgtmeantres >=r(p1) & wbgtmeantres<=r(p99)
replace wbgtmeantres = round(wbgtmeantres,0.1)
g ledbin = ledres>=0
g count=1
collapse (sum) count, by(wbgtmeantres led)
bys led: egen total = total(count)
g percentobs = count/total
twoway (scatter percentobs wbgtmeantres if led==0, mcolor(red) msymbol(circle_hollow) ) ///
(scatter percentobs wbgtmeantres if led==1,  mcolor(blue)) ///
, scheme(s2color) legend(order(1 "Without LED" 2 "With LED")) ///
xtitle("Wet Bulb Globe Temperature Residuals") ytitle("Percent of Observations", height(7)) ///
note("Residuals from regressions on budgeted efficiency and unit by year, unit by month, day of week, and line FE.", size(vsmall)) ///
caption("Temperature residuals are trimmed at the 1st and 99th. Temperature residual bins have width of .1", size(vsmall))
	graph export "tempdistbyled_wbgt.pdf", replace	
			graph export "LEDAppendixTableB.pdf", replace	

twoway (scatter percentobs wbgtmeantres if led==0, mcolor(black) msymbol(circle_hollow)) ///
(scatter percentobs wbgtmeantres if led==1,  mcolor(black)) ///
, scheme(s2color) legend(order(1 "Without LED" 2 "With LED")) ///
xtitle("Wet Bulb Globe Temperature Residuals") ytitle("Percent of Observations", height(7)) ///
note("Residuals from regressions on budgeted efficiency and unit by year, unit by month, day of week, and line FE.", size(vsmall)) ///
caption("Temperature residuals are trimmed at the 1st and 99th. Temperature residual bins have width of .1", size(vsmall))
	graph export "tempdistbyled_wbgt_bw.pdf", replace	

		graph export "LEDAppendixTableB_bw.pdf", replace	

restore
